\documentclass[10pt]{article}
\usepackage{amsmath}
\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex}

\newcommand{\problem}[2][No]{\noindent{\sc Problem #1.} #2}
\newcommand{\solution}{\noindent{\sc Solution.}}
\newcommand{\answer}[1]{\noindent{\sc Answer. }#1}
\newcommand{\complexity}{\noindent{\sc Complexity. }}
\newcommand{\reference}{\noindent{\sc Reference. }}

\newcommand{\given}{\,|\,}
\newcommand{\divides}{\,|\,}
\newcommand{\BigO}{\mathcal{O}}

\begin{document}

\problem[153]{Investigating Gaussian Integers.}

A \emph{Gaussian integer} is a complex number $a+bi$ such that both $a$ and $b$ are integers.

The regular integers are also Gaussian integers (with $b=0$). To distinguish them from Gaussian integers with $b \ne 0$, we call such integers \emph{rational integers}.

A Gaussian integer is called a divisor of a rational integer $n$ if the result is also a Gaussian integer. Note also that if the Gaussian Integer $(a+bi)$ is a divisor of a rational integer $n$, then its complex conjugate $(a-bi)$ is also a divisor of $n$.

The following is a table of all of the Gaussian integer divisors with positive real part for the first five positive rational integers:
\begin{center}
\begin{tabular}{c l c}
\hline
$n$ & Gaussian integer divisors & $s(n)$, Sum of divisors \\
\hline
1 & $1$ & 1 \\
2 & $1, 1+i, 1-i, 2$ & 5 \\
3 & $1, 3$ & 4 \\
4 & $1, 1+i, 1-i, 2, 2+2i, 2-2i, 4$ & 13 \\
5 & $1, 1+2i, 1-2i, 2+i, 2-i, 5$ & 12 \\
\hline
\end{tabular}
\end{center}

For divisors with positive real parts, we have: $\sum_{n=1}^5 s(n) = 35$.

It is also given that $\sum_{n=1}^{10^5} s(n)=17924657155$.

Find the value of $\sum_{n=1}^{10^8} s(n)$.

\solution

We split the problem into counting the real divisors and complex divisors separately. 

Let $R(N)$ be the sum of all real divisors of $n$ where $1 \le n \le N$. To find $R(N)$, we can enumerate all integers $1 \le k \le N$. For each $k$, there are $\lfloor N/k \rfloor$ multiples of $k$ below $N$. Each multiple of $k$ contributes $k$ to the sum. This gives an expression of $R(N)$ as below:
\[
R(N) = \sum_{k=1}^N k \left \lfloor \frac{N}{k} \right \rfloor .
\]
% Note the symmetry in the above expression: for every $k>\lfloor \sqrt{N} \rfloor$, $\lfloor N/k \rfloor < \lfloor \sqrt{N} \rfloor$. That is, the products come in pairs. We can therefore rewrite $R(N)$ as
% \[
% R(N) = \lfloor \sqrt{N} \rfloor^2 + 2 \left( \sum_{k=1}^{\sqrt{N}-1} k \left \lfloor \frac{N}{k} \right \rfloor \right).
% \]
% This expression is computationally much simpler than the original one.
Note, however, that every pair of divisors $(b,a)$ where $b>a$ can be mapped to $(a,b)$ where $a<b$. Therefore, we only need to enumerate pairs $(a,b)$ where $a \le b$. In this approach, we enumerate $a$ from 1 to $\lfloor \sqrt{N} \rfloor$. For each $a$, we enumerate $b$ from $a$ to $\lfloor N/a \rfloor$. Each pair of $(a,b)$ generates a product $ab = n \le N$. This provides another expression of the sum of real divisors as
\begin{align}
R(N) 
&= \sum_{a=1}^{\lfloor \sqrt{N} \rfloor} \left[a + \sum_{b=a+1}^{\lfloor N/a \rfloor} (a+b) \right] \notag \\
&= \sum_{a=1}^{\lfloor \sqrt{N} \rfloor} \left[ \frac{(3a + \lfloor N/a \rfloor + 1)(\lfloor N/a \rfloor - a)}{2} + a \right] \notag
\end{align}
This expression is computationally much simpler than the original one.

To find the sum of all complex divisors of each $n$ where $1 \le n \le N$, we enumerate all complex pairs $(a+bi, c-di)$ whose product $(ac+bd) + (bc-ad)i$ is a real integer less than or equal to $N$. This leads to two conditions:
\begin{equation}
ac+bd \le N \label{eq:one}
\end{equation}
and
\begin{equation}
bc-ad = 0 \label{eq:two}
\end{equation}
Note that when we find a divisor $(a+bi)$, we also find three other divisors $(a-bi)$, $(b+ai)$, $(b-ai)$, because
\begin{align}
(a-bi)(c+di) &= (ac+bd) + (ad-bc)i = n \notag \\
(b+ai)(d-ci) &= (ac+bd) + (ad-bc)i = n \notag \\
(b-ai)(d+ci) &= (ac+bd) + (bc-ad)i = n .\notag
\end{align}
Therefore we only need to enumerate divisors of the form $a+bi$ where $a \le b$; any other divisor can be mapped to a divisor of this form. If $a<b$, the group of four divisors add to $2(a+b)$. If $a=b$, the pair of two divisors add to $(a+b)$.

Rearrange equation \eqref{eq:two} as
\[
\frac{a}{b} = \frac{c}{d} = \frac{p}{q}
\]
where $p/q$ is the reduced form of the fraction. Let $u = \gcd(a,b)$ and $v = \gcd(c,d)$, and we can write the fractions as
\[
\frac{a}{b} = \frac{pu}{qu} \text{ and } \frac{c}{d} = \frac{pv}{qv}
\]
where $u,v \ge 1$. Substitute this into equation \eqref{eq:one}, we get 
\[
ac+bd = (p^2+q^2)uv \le N .
\]

The algorithm is as follows. Given coprime integers $p<q$, we enumerate $a+bi = u(p+qi)$ for $u$ from 1 to $\lfloor N/(p^2+q^2) \rfloor$. For each $u$, there are $\lfloor N/[(p^2+q^2)u] \rfloor$ products that are less than or equal to $N$. Each product contributes $2(a+b)$ to the sum of divisors. Therefore, the sum of divisors generated from coprime pairs $p<q$ is
\begin{align}
C_{p,q}(N) 
&= \sum_{u=1}^{\lfloor N/(p^2+q^2) \rfloor} \left\lfloor \frac{N}{u(p^2+q^2)} \right\rfloor 2u(p+q) \notag \\
&= 2(p+q) \sum_{u=1}^{\lfloor N/(p^2+q^2) \rfloor} \left\lfloor \frac{\lfloor N/(p^2+q^2) \rfloor}{u} \right\rfloor u \notag \\
&= 2(p+q) R(\lfloor N/(p^2+q^2) \rfloor) .\notag
\end{align}
The only coprime pair where $p=q$ is $(1,1)$. The sum generated from this pair is
\[
C_{1,1}(N) 
= \sum_{u=1}^{\lfloor N/2 \rfloor} \left\lfloor \frac{N}{2u} \right\rfloor 2u
= 2 R( \lfloor N/2 \rfloor ).
\]

Finally, we need a way to generate coprime integer pairs $(p,q)$ where $p<q$ and $p^2+q^2 \le N$. If we interpret $(p,q)$ as the coordinates of a point on the plane, then this is equivalent to enumerating all integer points within the circle $x^2+y^2 = N$.

We make use of the properties of the \emph{farey sequence} of order $k$, denoted $F_k$, to do the enumeration. The steps are:

Step 1. Start with $k=\lfloor \sqrt{N-1} \rfloor$, and let $p/q=1/k$. This is the smallest term of the farey sequence $F_{k}$.

Step 2. If $p^2+q^2 \le N$, compute the sum generated from this coprime pair. Then let $p/q$ be the next term in the farey sequence.

Step 3. If $p^2+q^2 > N$, this means all the remaining fractions with denominator $q$ will also be outside the circle, so we don't need to enumerate other fractions with this denominator any more. Therefore, we decrease the order $k$ by one.

Step 4. If $q=1$, then finish. Otherwise, go to Step 2.

Note 1. A technical detail here is how to compute the next term in a farey sequence, given the previous two terms. It can be shown that if the previous terms are $p'/q'$ and $p/q$, then the next term is 
\[
\frac{\lfloor (r + q')/q \rfloor p - p'}{\lfloor (r + q')/q \rfloor q - q'}
\]

Note 2. An intensive operation in the above algorithm is computing 
\[
R(\lfloor N/(p^2+q^2) \rfloor)
\]
for all coprime pairs $(p,q)$. It turns out that for large $(p,q)$, $\lfloor N/(p^2+q^2) \rfloor$ are small and often have the same value. We can therefore improve the running time (at the cost of more memory usage) by caching $R(n)$ for small $n$. Here we choose to cache $1 \le n \le \lfloor \sqrt{N} \rfloor$.

% From Wikipedia, we find the following generating formula:

% Starting with $(m,n) = (2,1)$ or $(3,1)$. For each step, three branches are generated:

% Branch 1: $m' = 2m - n, n' = m$

% Branch 2: $m' = 2m + n, n' = m$

% Branch 3: $m' = m + 2n, n' = n$

% The process is continued recursively. Each generated pair $(m,n)$ satisfies $m > n$. Finally we add the pair $(1,1)$.

\answer{17971254122360635}

\complexity

First we estimate the number of coprime pairs $(p,q)$ such that $p^2+q^2<n$ as a function of $n$. The density of coprime pairs is $6/\pi^2$. The area of a quarter circle of radius $\sqrt{n}$ is $\pi n/4$. So the number of coprime pairs in a quarter circle of radius $\sqrt{n}$ is about $3 n/(2\pi) \approx 0.48 n$, which is of the order $\BigO(n)$.

In the non-cached version of the algorithm, the space complexity is $\BigO(1)$. In terms of time complexity, computing $R(N)$ costs $\BigO(\sqrt{N})$. Computing $C_{p,q}(N)$ for a given pair of $(p,q)$ costs $\BigO(\sqrt{N/(p^2+q^2)})$ which is no more than $\BigO(\sqrt{N})$. Since the number of coprime pairs is $\BigO(n)$, the time complexity is $\BigO(n^{1.5})$.

In the cached version of the algorithm, where we cache the value of $R(n)$ for $1 \le n \le \sqrt{N}$, the space complexity is obviously $\BigO(\sqrt{n})$. In terms of time complexity, it costs $\BigO(n^{0.75})$ to build the cache. Then, given coprime pair $(p,q)$, if $p^2+q^2 \ge \sqrt{N}$, then $\lfloor N/(p^2+q^2) \rfloor \le \sqrt{N}$ and the value of $R(\cdot)$ can be retrieved from the cache and costs $\BigO(1)$. There are $\BigO(N)$ such pairs. If $p^2+q^2 < \sqrt{N}$, it costs $\BigO(\sqrt{N})$ to compute $R(\lfloor N/(p^2+q^2) \rfloor)$. There are $\BigO(n^{0.25})$ such pairs. Overall, the time complexity is $\BigO(N)$.

The complexity is summarized in the following table:
\begin{center}
\begin{tabular}{l | c c}
\hline 
Algorithm  & Time complexity  & Space complexity \\ 
\hline
Non-cached & $\BigO(N^{3/2})$ & $\BigO(1)$ \\
Cached     & $\BigO(N)$       & $\BigO(N^{1/2})$ \\
\hline
\end{tabular}
\end{center}

\reference

http://en.wikipedia.org/wiki/Coprime

http://en.wikipedia.org/wiki/Farey\_sequence

\end{document} 